function dz = diff_Reg1(t,z)

% Differential equations in Regime 1, where only fossil fuel is  used. iR>>0, n>0.

global par

k = z(1);
lambda = z(2);
f = z(3);
kR = z(4);
S = z(5);
N = z(6);
sigma = z(7);


dz = zeros(7,1);

%%% calculate the variables 

F = par.F0-par.a*f;
denom_pe = par.a*k + F;
pe = lambda * (par.A - par.delta) / denom_pe;

g = par.alpha0 + par.alpha1/(par.Sbar-par.alpha2/(par.alpha3+N)-S);
gdS = alpha1*(alpha3+N)^2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gdN = -alpha1*alpha2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gd2N = 2*alpha1*alpha2*(Sbar-S)/((Sbar-S)*(alpha3+N)-alpha2)^3;
gdSN = -2*alpha1*alpha2*(alpha3+N)/((Sbar-S)*(alpha3+N)-alpha2)^3;

           
c = lambda^(-1/par.gamma);

Q = par.Q0.*exp(par.popgr*t);

%%% change some denotation to simplify the expression
combo_qR = par.muR*(sigma*Q-lambda*g);
                             


    
    LHS = zeros(4,4);
    RHS = zeros(4,1);  
         
    LHS(1,:) = ones(1,4);
    LHS(2,:) = [(par.F0 - par.a*f), -par.a*k, -par.muR, 0];
   
    LHS(3,:) = [combo_qR+par.A*lambda, -(combo_qR+par.delta*lambda), ...
               0, -(par.A-par.delta)*k*lambda/kR];
    LHS(4,:) = [0,-2*F*gdN^2, par.muR*gdN^2, (par.A-par.delta)*gd2N];
    
    RHS(1) = par.A*k - g*par.muR*kR - c;
    RHS(2) = 0;
    
    RHS(3) = par.delta*k*(combo_qR+par.A*lambda)+ (par.A-par.delta)* ...
             sigma*Q*k*(par.muR-par.popgr*(kR*gdN)^(-1));
    RHS(4) = par.delta*F*gdN^2+(par.A-par.delta)*gdSN*par.muR*kR*Q;
    
    
    
    sol = RHS\LHS;
    
    inv = sol(1);
    iF = sol(2);
    iR = sol(3);
    n = sol(4);
   
    

dz(1) = inv - par.delta*k;
dz(2) = (par.beta + par.delta - par.A) * lambda + F * pe;
dz(3) = iF;
dz(4) = iR-par.delta*kR;
dz(5) = Q*par.muR*kR;
dz(6) = n;
dz(7) = par.beta*sigma + lambda*par.muR*kR*gdS;
